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ABSTRACT 


The governing partial differential equations of a 
classical shell theory are reduced to a set of eight first 
order ordinary differential equations. A forward numerical 
integration process is used to obtain influence coefficients and 
particular solutions for shells of revolution of general geometric 


configuration subjected to arbitrary types of loading. 


Standard stiffness methods of structural analysis 
are employed to obtain displacements and stress resultants 
everywhere within a complex shell structure. A computer program 
is developed to perform the analysis based on the theory presented. 
Example problems are selected to test the accuracy of the method. 
Excellent results, when compared with known solutions, are 


achieved. 
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circumferential and meridional twisting moments per 

unit length, respectively. 

harmonic number 

normal in-plane forces per meridional and circumferential 
unit length, respectively. 

in-plane shear forces per meridional and circumferential 
unit length, respectively. 

intensity of load components in the direction s, z and 

8 respectively. 

transverse shear forces per meridional and circumferential 


unit length, respectively. 
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column vector arises from integrating the inhomogeneous 
terms in the eight fundamental equations (Eq. 3.14). 
radius of curvature of parallel circles. 

radius of curvature of meridian. 

length of the normal between any point on the middle 
surface and the axis of revolution. 

derivative of r,; with respect to the coordinate s. 
curvature of paraltiel circles: 

fivst principle curvature = 1/71. 

second principle curvature = 1/rp. 

coordinate measures the distance along the meridian 

of the shell. 

effective transverse shearing force per circumferential 
unit length defined in Eq. 2.6.1. 

thickness of the shell. 

change in temperature from arbitrary level. 

effective tangential shearing force per circumferential 
unit .Jéeneth defined in Eq. 2.6.2. 

temperature terms defined in Eq. A.19.1. 

temperature terms defined in Eq. A.19.2. 

change in temperature from arbitrary level x 

measured at the exterior and interior face of the shell, 
respectively. 

displacement component in the circumferential direction. 


displacement component in the meridional direction. 
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displacement component in the rotation direction. 
displacement components of a point at distance z 
from the middle surface. 

coordinate measures the distance along the axis of 
revolution. 

column vector of the eight fundamental variables 
CEd.t 54.0). 

derivative of Ly page with respect to the coordinate s. 
coefficient of thermal expansion. 

rotation of the meridian due to deformation. 

specific weight of the shell material. 

shear strain. 

meridional strain. 

hoop strain. 

coordinate measures the angle in the circumferential 
direction. 

Poisson's ratio. 

normal stresses in the meridional and circumferential 
direction, respectively. 

shearing stresses. 

coordinate measures the angle between any point on 


the middle surface and the axis of revolution: 
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CHAPTER 1 


INTRODUCTION 


i haa Introductory Remarks 


Shells of revolution are important structural elements. 
Many structures such as storage tanks, pressure vessels, silos, 
chimneys and towers are composed of either a single shell unit or 
an assemblage consisting of different types of shells. Their 
behaviour allows the shell thickness to be reduced to a minimum 
and their advantageous shapes permit more modern architectural 


CORCEDES|. 


Although the governing field equations have been known 
for many years, cases where analytical solutions can be obtained 
are relatively scarce and are restricted to simple forms of 
geometry, boundaries and loads. The determination of the forces 
and deformations in shells constitutes a difficult problem in 
the theory of elasticity, owing to the complexity of the mathe- 
matical equations involved. For conditions in which the analy- 
tical solution is complex, or is unknown, the application of 
numerical methods with the aid of a digital computer has proven 
to be useful and efficient. This approach allows the solution 


for generalized geometric configurations and loadings of shells 


of revolution. 
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1.2 Purpose of the Study 


The objectives of this thesis are: 

1) to develop a technique for evaluating the stiffness 
influence coefficients of any arbitrary element of a 
shell of revolution of general geometric configuration 
by means of a direct numerical integration method. 

2) to employ the standard stiffness method of structural 
analysis to analyze assemblages of such elements. 

3) to demonstrate the capability of the method to treat 
arbitrary surface loadings and thermal gradients, 


including line and concentrated loads. 


1.3 Types of Shells of Revolution 


The position of a point in a shell of revolution can 
be given by three orthogonal coordinates s, 9 and z (See Appendix 
A, Sect. A.1, for definitions). The shape of the shell is 
determined by specifying the two principal radii of curvature 
Yr} and ro of the middle surface and the thickness of the shell 
(Fig. A.1). Instead of r, it is sometimes convenient to use the 
distance r from a point on the middle surface to the axis of 


revolution (Fig. A.1) where 


r = r,sind i Ba | 


in which > is the angular distance of the point under consideration 


from the axis of revolution. The generating curve of the middle 
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surface is defined by the equation 


r= r(x) 


where r(x) represent the radius r as a function of the distance 
measured along the axis of revolution, x. Therefore the 
principal radii of curvature can be determined by the following 


two equations 
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The general shape of any type of shell of revolution 
(Fig. 1.1 towh6) Is ‘characterized by particular worms Om Eqis. 


as follows: 


aieeior plates (hig. 1.1) Cie oe 
ro ee 
d5=) 0 
Bb) for spheres (Fig. 1.2) ty <8 
fe =e 
Oe tor cylinders, (Pig.yi.3) mr, = © 
Ee oe 
od = m/2 
d) for cones (Fig. 1.4) r, = © 
lf ere 3 
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S) for coroids® (Fig .© 145) OY pa Epa. 
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roe = ER 1s 


f) for ellipsoids (minor axis coincides with the 


axis of revolution) (Fig. 1.6) 
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h) for arbitrary shell elements for which the form of the 
middle surface cannot be expressed as a closed form 
function (Pig: 1.3), one can describe r,, 4, and 1 at 
discrete points along the meridian and interpolate 


numerically. 


1.4 Loadings 


Applied surface loads at any point of the shell can be 
resolved into three components in the three orthogonal directions 
sg, 9 and zg. This load may vary in the direction along the meridian 
as well as in the circumferential direction of the shell. 

Therefore a load component may be written as a function of the 


coordinates s and 6 in the following form 


= 


P, = Fi,(s,6) 1.5. 
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where P. is the magnitude of the load at the point under considera- 
tion in the direction i(i = s,8,z) and Fis is the function 
representing the applied load. For the special case of axisymmetric 


loading P. is independent of 9, thus 
fg ao Eos 


For non-axisymmetrical loading, the classical method 
of analysis is to expand this load into a Fourier series, 
analyze for each harmonic separately and superimpose the effects 
[10,27]. The number of terms considered in this series must be 


sufficient to give the desired degree of convergence. 


Thermal loading can be considered in either case by 
algebraically adding the strain due to thermal expansion to the 


strain due to the surface loading in the stress-strain equations. 


1.5 Shell Theory 


Shell theories of various degrees of complexity may be 
derived, depending upon the degree to which the theory of linear 


elasticity is simplified. 


In all cases one begins by reducing the three-dimensional 
shell problem to a two-dimensional problem expressed in terms of 
the deformation of the middle surface of the shell. Further 
simplifications establish various shell theories which may be 
classified into different categories [3]. Such categories are 


based on the terms that are retained in the strain and stress- 
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displacement equations with respect to the thickness 


coordinate. 


The second order approximation theory for shells of 


revolution was presented in 1932 by Flugge [10] and based upon 


the following assumptions: 


1) 


2) 


3) 


Consistent with the formulation of the classical theory 
of elasticity, strains and displacements under loads 
are small enough so that changes in geometry of the 
shell will not alter the equations of static equili- 
brium of the shell (i.e., equations of equilibrium are 
written in the undeformed configuration). 

The components of stress normal to the middle surface 
are small compared to the other components of stress 
and may be neglected in the stress-strain relationships 
(i.e., the material may be considered to be in a plane 
stress condition). 

Points on lines normal to the middle surface before 
deformation remain on straight lines normal to the 
middle surface after deformations (i.e., deformations 


of the shell due to the radial shears are neglected). 


Contrary to other theories, Flugge's theory did not entirely 


neglect the ratio of the thickness to the radius of curvature 


(except for an occasional dropping of the fifth and higher order 


terms) in the stress resultant equations and in the strain- 


displacement relationships [10, pg. 320]. 
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Applications of this theory have generally been 
restricted to circular cylindrical shapes, for which some 
solutions have been obtained [10,16,19] and are considered as 
standards against which other simplified theories may be compared 


[oe £9)". 


The fundamental equations of Flugge's theory, upon 


which this study is based, are presented in Appendix A. 


1.6 Methods of Analysis 


A structure usually consists of an assemblage of many 
parts and tends to be complex in nature. Generally, the true 
structure must be replaced by an idealized approximation, or 
model, suitable for mathematical analysis. Structural analysis 
for shells may conveniently be carried out by matrix methods using 


influence coefficients. 


In the literature, the analysis of symmetrically loaded 
shells of revolution is classically performed using flexibility 
influence coefficients [3,4,21]. These have been given, in an 
explicit manner, for very limited number of types of shells of 
revolution, such as, cylinders, spheres and cones, of uniform 


thickness. 


For symmetrically loaded shells, the membrane solution, 
which represents the momentless state of stress in the shell, 
may approximate the particular solution which satisfies the 


general differential equation of the shell. Using the membrane 
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solution, a flexibility method of analysis can be performed by 
satisfying the continuity requirements at the joints at which 
the elements are connected. When the forces and displacements 
at these joints are known, the conditions within each element 
may be determined from the solution of the differential equation 


of the element. 


For the case of arbitrary shells of revolutions, under 
arbitrary systems of loadings, the corresponding expressions for 
influence coefficients are unknowns. The analytical method of 
obtaining these expressions would involve the solution of eighth- 
order differential equations expressed in terms of the geometry 
of the shell surface and the physical constants. This method is 
difficult and complicated, even for a simple geometry such as a 
cylinder, and it is generally impossible for the case of an 


arbitrary shell under arbitrary loads. However, the latter 


case is the rule rather than the exception in modern architecture. 


As a consequence of the availability of electronic 
computers, and the increasing familiarity of engineers with this 
computational tool, the application of numerical methods of 


analysis to shell problems has become more attractive. 


Two numerical methods for the analysis of shells of 
revolution with arbitrary configurations have received extensive 


treatment in the technical literature. The first method is the 


finite difference method, which consists essentially of the direct 


replacement of the derivatives which appear in the governing 
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differential equations by finite difference approximations. 

This method transforms the differential equations into a system 
of algebraic equations which may be solved by an iterative 
procedure [26] or by means of matrix methods [6]. The method is 
quite general in application. Replacement of the derivatives by 
finite difference approximations may be undertaken at any level 
in the basic formulation of the shell problem. However, it is 
difficult to introduce boundary conditions into the problem. It 
also becomes cumbersome when attempting to satisfy equations 


involving high order derivatives [9,12]. 


The second method is the finite element method. In 
this method the displacement of each element into which a shell 
is subdivided are represented by an approximation [9,11,13,20,24]. 
The most common practice is to represent each shell by a series 
of short conical shell elements of uniform thickness. The 
variation in thickness along the generator of the shell, can be 
accounted for by considering the average thickness of each 
element [20]. Once these short segments have been defined, the 
problem becomes one of analyzing a shell that is an assembly of 
many short conical shells. Stiffness influence coefficients are 
evaluated using energy methods. Conditions of continuity are 
then applied at the boundaries of each segment to evaluate the 


forces and deformations. 


An alternative method of numerical analysis based on 


direct integration of the shell equations, as proposed by 


2 
] 


: f Nop } . we 
- ; ; i» ; i ' 
(wo : : ui oa 
tt 7 7 i Gs or ane i 
oa) On } t y ‘a i 
oe ts, i” | ' t ; 

[hie "al 7 

a, 7 14 A 

7 1 i . 


ie a, eA okt sciobipndsorga Sanabaatiind BS 


iF 
m meleere 5 want” stata slips: cniahionan oy sven _— 
| ner er ‘ie we | te) fea) ae saa nok taupe ‘4 s d | 
Se hgiecunuees Fea whorl 2h site abisive ee; 30 pei 
ae 4 Sercy ee 10 (apne: thei 3 Rvciamnaa’ ait 7 
| oT rie nsdecy sii ad “go bashiamornae 
ed 4 EF) Liwpeiee  vmies oe Pied ‘fi ‘Somssrus a 
+) . sldete saa © ‘ob 3 Y apiet se? i wai" V7 aia lpi Yar 
| Maat js) yite ba: 63 ou £iqpg 33k ‘a dows : it 
\ feb ie] 
Gi) )-p0na joa! orpmets = dnp ha 
Lede .S Kobe eet, dnshs is dome ms derauate alee oa 
; AP OEE L100) idawectaue tes ta Se biped: 
BILE, role L tees Hae tee oe aa sm 
| | siiiin eee Na “ei to Hobmmisetesi | j 
uh : ie asks) dttoie oti 14) 7 .98Sona ld, BeuTh eee, 1 
rs | | i309 QB; iid & - tis oily ory Deana ys 70 
: eles sbedu iat Root ath, fn Se nats sti ata Ah 
. ahs i te vilowzan ne ak tee aes a eh! Sct vite eo. ponooga 


as. Sins ste, Cae ») pels ee hela ig Loma Thais ane . ce 


‘ i i) eriyanaee ie Ser Lbady: al tn ei eth be: 


vate edits: et a alle fans er ietimgad e . 
nar 4 ‘ Pe oly, a we bettas 


me ho ee ell _emnatseminiies bee 
ie 0 faa acoilend shaban. ty Ss Pasenee 
Ly om dasa 26 dienes Liste oe te | 


Ae he | ae. 
: ak aa : — i i : ie 3 


i % 
rv ’ ® ma © 


LO; 


Goldberg, et al [12] for spherical domes and by Iyer and 
Simmonds [17] for conical shells, is presented for general 


shells of revolution in this thesis. 


In the application of this method to the analysis 
of shells of revolution under arbitrary loadings, the governing 
partial differential equations of a consistent shell theory are 
expanded by means of Fourier series. Differentiation with 
respect to the colatitude coordinate 0 may be performed to 
transform the governing equations to ordinary differential 
equations. These equations are reduced to a set of eight first- 
order ordinary differential equations involving eight intrinsic 
dependent variables and their derivatives. The intrinsic variables 
are the three components of displacement at the middle surface, 
the rotation of the tangent to the meridian, the membrane normal 
force in the meridian direction, the moment acting on circum- 
ferential sections and two modified shear terms in the directions 


perpendicular and tangent to the meridian. 


This system of equations can be integrated in a stepwise 
fashion across a given element. By performing matrix operations, 
which will be described in Chapter 3, the stiffness influence 
coefficients can be calculated. Fixed end forces can be obtained 
and stiffness analysis may then be performed to determine the 
conditions at the element boundaries. These conditions can be 
used to determine the forces and deformation everywhere within 


the shell element. 
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i.7 “Outline of Contents 


The governing equations of the classical theory of 
shells of revolution are reduced to a set of eight basic equations 
in Chapter 2. In Chapter 3, the solution technique for the basic 
equations is described. Standard stiffness methods for segmented 
shell structures are outlined. Two example problems are handled 
in Chapter 4. General types of loadings are considered. In 
Chapter 5, limitations of the technique are discussed. Conclusions 
are drawn and possible future developments are outlined in 


Chapter 6. 
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Fig. 1.8 
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CHAPTER 2 


BASIC EQUATIONS 


2.8) Introduction 


In this chapter the governing field equations of a 
shell of revolution, as derived in Appendix A, are expanded using 
Fourier series. Modified shear terms are introduced into the 
governing equations to eliminate the inplane shearing force in 
the circumferential direction and the meridional twisting 
moment, which appear at the boundaries s = constant. By performing 
matrix operations to eliminate the stress resultants in the 
circumferential direction (secondary stress resultants) the 
governing field equations are reduced to a set of eight first 
order equations involving four stress resultants (primary stress 


resultants) and four displacements and their derivatives. 


2,2 fourier Series 


The load components and temperature, being arbitrary 


functions of s and 8, may always be represented in the form 


P =  P (s)cosn6+ £ P_ (s) sinné Ziel, 
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(oe) 
T = X T_ (s)cosn0 + £ T (s) sinné 
n n 
n=0 n=l 
where Py Py and P. are defined in Appendix A, T represents 
temperature, and the functions of s on the right hand side are 


Fourier coefficients [10,27]. The corresponding stress 


resultants and displacements may be expressed as 


N = > N_(s)cosn8 + XN (s) sinné 
Ss sn sn 
n=0 n=1 
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n=0 n=l 
Nog = bx Ma (s) cosnO + : N oon (s) sin nO 
n=0 n=1 
No, = » ee (s) cosnO + : Noon (s) sin n9 
n=0 n=1 
Qg = D <a (s) cosn0 + 2% Cie (s) sin nd 
n=0 n=l 
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Q% = : Qon (s) cosnO0+ 2% Qon (s) sin nd 
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n=0 n=1 
= M 6+ XM (s) sin nO 
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where the variables on the left hand side are defined in Appendix A, 
Sects. A.2 and A.3, and where the term (s) indicates that the 
variable coefficients with subscript n are functions of the 
coordinate s only and n is the harmonic number. The first and 

the second series in each expression represent the portions of 

the variables which are, respectively, symmetric and anti- 
symmetric with respect to the meridian passing through the line 


0=0. 


For an arbitrary applied load expressed as a Fourier 
Sel ros worm ciel Order iN (EGS. 2 a *to 2.9.3)" there are 2Ner 
terms that’ represent each component of ‘load; (n = 0, 15° 25...,N) 
tor the symmetric series and (n= 1} 2, <%., N) for ‘the anti- 
symmetric series. For each value of n the (s dependent) 
variables with the subscript n from each series can be entered 
in the governing equations of shells of revolution of Appendix A, 
because the sequences cosn9 and sinnO are linearly independent. 
Differentiations with respect to 9 can be performed and the 


terms grouped according to the common factors, cosn9 or sinnQ. 


Since the coefficient of each of these factors must be zero, each 


x j ~ 7 Aon | an 

re ' to eR oo 
it ; ‘ OF ne 
"i H _ . a ; 
q ; ] iin Marca ty: 
a * bel ye 1 is ie x iJ - i 
, Di 'Si% ‘gie te) G wae” Tis) 

ie , 


te | Shinke Coy ‘- e + anion 1) - ra r a 7 # 


¥ 
~ py 1 fi 


Pawel : ye . ae > 


4 ft 
(COR fa ; rt 


Rises Ooi inks (ays ow ish han Ronete) Saha aes 
le : L=f; ry oe ie 


‘yA. Bibasgca: ol, beri $4 oTh Lite bitsy si me wd 

a ae ra a iat Bee Be: Hy SHS be Cayo miat ‘eit qe 151 Na ‘ba ser bes: wna 
. : } | ay 46 3 neds air er bow 74 G4 avon tie agi 3 irate eee " ald 
a | Wak g245% ed iedmae. Shore 1 orld Ba? a ike! wf, Le | 
a eres et ics : abi pte wits POE e “ sake ia | oy ‘bite | 


re taf ast bs iy | he} reveenyy = . ais ie: ¥ t + Haul eat" al Bb ai : ys : 


ye 


pa) 


an 


si neh by eek? od deooids 4nteasg nathiaom aad ci? 9g darcy Hahw. 


; ' i i 7 ; P > wa! . iu 


nee RStano OSes eseaatqk> bibs So bts spect ma i, 


fous sts) otoi1", (fabs parteiak : apo), W abe, al 


eases tT) 0) <\ a) Pb fl AO “ptogey > co dips 
Mebshe hdd TR Ce a a Ges at aiktus) wb se 
| . a ae Rise a) Oh Buds paidawel die gt i ested b 


bezels ‘Sd 62 es i138 | ae! gs ” 4: bisatue) sik peat nae 


a A, stbooaih They) mopaulaead bc abides ty ped sepa, 


o Bi Adina. aicath Be iy ie ba Ou as sa aomne phere 
7 | : | 


“ttt on thorgetns a0 aso us ‘ed losiggee thie '¢ 


im, 


vies ni to. yt Zon pete oa | phates = ay, pean 
otis aa Feuim b20a0h% Seti, 4G; feos bad aes be 


| a 
pure 4 mt ) se , A y os 


ft 


Ved 


factor produces a separate equation. For example, Eqn. A.5.1 is 


rN + Ne Aa Bp a 
fi cos > is Qs cos ON, = OF + rP 0 Z asic 


for which, for any n, the cosine terms become 


rN’ _cosn0 + cosdN__cosnO + nNg __ cos nO 
sn sn sn 


r 
- cosON, cos no - r Qon COS n6 + rP cos nO = O Dui dhs 


which, upon factoring out the common term, yields 


i 
- cos Ny =a) 


rN’ COS ON a. tae N 
sn sn Le Sst 
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+ rP = 0 2° 


If the Fourier expansions of Eqs. 2.1 and 2.2 are 
truncated such that 0 < n < N, the governing system of equations 
is replaced by (2N+1) ordinary differential systems of equations, 
each equation of the type illustrated by Eq. 2.3.3, and an 
analysis is carried out to obtain the solutions for the s 
dependent coefficients corresponding to each value of n. In 
the general case there are twenty-six variable coefficients 


associated with each n > 0 in Eqns. 2.2, thirteen associated 
with each of the factors. 

Finally, since linear stress-displacement relations 
were assumed, the principal of superposition is valid and one can 


superimpose the QN+1) solutions so obtained. 
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The nt) set of equations can be written as follows. 


The six equilibrium equations are obtained from Eqs. A.5 as 
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where Rj and Ro are the principal curvatures of the shell, 


defined as the reciprocals of the radii of curvature r, and r,, 


respectively, 


The 


from Eqs. Al/7 
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and R is the curvature of the parallel circle. 
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eight stress-displacement equations are obtained 


as 
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derivatives of the displacement component W. 
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the thirteen unknown functions (the variable Fourier series 
coefficients). As such the system of equations is overspecified. 
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remaining thirteen equations allow the solution for the thirteen 
functions associated with each trignometric functions of the 


Fourier series expansion. 


22 Natural Boundary Conditions 


According. to the classical theory of shells, the 
quantities which appear in the natural boundary conditions on an 
edge s = constant of a shell of revolution are the four displace- 
ments 8, W, V, U and the corresponding four forces M.> S No» 
T.: 

The forces pes T. are the transverse and tangential 
effective shears which are commonly known as Kirchhoff's shears [27]. 
The Kirchhoff shears are work-equivalent forces associated with 
the displacements W and U, and these effective shears replace the 
stress resultants ee Nog and Mog: Such replacement is essential 
to describe simple boundary conditions [10,27]. Although the 
expressions for these forces may be derived rigorously from 
variational principles, the following physical interpretation, 
due to Kelvin and Tait [27}; pe:45li, is more instructive. If 
the twisting moment Mog acting on an infinitesimally small 
element of the shell is replaced by a statically equivalent 


force as shown in Fig. 2.1, one can write fromstatics the 


expression for the effective transverse shearing force as 1274 
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and the effective tangential shearing force as 
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Writing Eqs. 2.6.1 and 2.6.2 in Fourier series terms, and 


separating the linearly independent terms, we may define 
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Using the geometrical relations in Eqs. A.2, the derivatives 
of these forces with respect to the coordinate s may be written 
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2.4 Reduction of the Governing Equations 


It can be seen that the stress resultants Qg » Nos M, 


M and N do not enter into any boundary conditions on an 
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where the form of the functions F, to Fy, is obtained from Eqs. 
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symbolic form, 
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of curvature (i.e., ri = ro), the sixth equation of equilibrium 
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eliminate them in terms of the other stress resultants and 


evaluate them after the solution of the governing equations. 


For convenience, the subscript n in Fourier coefficients 
will be omitted and the governing equations (Eqs. 2.4 and 2.5) 
can be written, for each set of equations, in the following 
symbolic form. The equilibrium equations (Eqs. 2.4) may be 


written symbolically as: 
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The stress resultant-displacement equations (Eqs. 2.5) may be 


written symbolically as: 
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Appendix A, Sect. A.5 and as noted in Sect. 2.3. 


In solving the above equations, the displacements 
which may be physically imposed on the boundary of a shell are 
8, W, V and U. The external forces associated with these 
displacements (in a work-equivalent sense) are M.; ape N, and 
T.: The first set of variables are specified in the case of 
displacement boundary conditions while the latter set is specified 
in the case of mechanical or force boundary conditions. In general 
a combination of four of these eight quantities must be specified 
to properly define a boundary condition, provided that if one of 
the displacements is specified the associated force should not 
be specified and vice-versa. Since the remaining five of the 
thirteen variables cannot be determined at the boundary, it is 
desirable to eliminate them from the governing equations. The 
objective in the following is, therefore, to reduce the thirteen 
governing equations (Eqs. 2.9.1 to 249.5 and Eqs: 2.10)-to a set 
of eight governing equations in terms of the eight variables 
which may arise in the specification of the boundary conditions. 
It is also desirable to keep the order of these differential 
equations to a minimum in order to facilitate the numerical 


solution. The final set of equations will be of first order. 
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Heingekas. 2.9.2 and 2.9.5 to eliminate Q% and then substituting 
for Mo5> Nog? N’ <6 by means of Eqs. Z27NI3 and 2.11.5, .one can 
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Eqs. 2.12, which represent the equilibrium equations, may now 
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5° = F 
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where the loading terms in these equations are separated in the 
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The coefficients of the matrices [Bl], [B2] and {B3} are defined 
any Table. 2:1 

In Eq. 2.14, the stress resultants have been separated 
into the vector oes whose components are desired in the final 
formulation, and the vector <Fo? > which remains to be eliminated 
from the formulation. 

Let us now turn our attention to the displacement 
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B°, V°, can be converted to the following two equations 
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Finally, Eq. 2.10.3, after eliminating N. by meansoor vq. 2.11.8 
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{p°} = [Al a2] f: + {a3} 2.18 
SS 
where 
<D> = < §6 W V LU Fog ae IE 
SPoa = <P “6s W° As ir, Ck are 


The coefficients of the matrices [Al], [A2] and the column 


vector {A3} are shown in Tables 2.2. 


The stress resultants appearing in the vector ig of 
Eq. 2.15.3 can be expressed in terms of the fundamental displace- 
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Therefore, for the applied loads, stresses and displacements 

being periodic functions of the coordinate 6, the general governing 
partial differential equations, Eqs. A.5 and Eqs. A.17, have been 
transformed into a system of eight first order ordinary differential 
equations for each harmonic member n. These equations relate, 

at any point, the eight fundamental dependent variables, that 

appear in the natural boundary conditions of shells of revolution, 


and their derivatives with respect to the independent variable s. 


As was seen, the reduction of the shell equation into 
this form involves only straight forward algebraic manipulations. 
The equations have been put in this form to facilitate numerical 
integration. They can be integrated with the aid of an appropriate 
integration technique, such as the Runge-Kutta process. Since 
a digital computer is needed to perform this integration, the 
reduction of the equations to the final form, as in Eq. 2.24, 
is not necessary and can be performed by feeding the three sets 


of equations (Eqs. 2.14, 2.18 and 2.20) into the computer. 


It can be seen that the elements of the matrix GA in Eq. 
2.24 are dependent only on the shell thickness, the physical 
constants and the coordinate s. It should be noted that in the 
case of shell of revolution under symmetrical loading conditions 
(i.e., n = 0) the system of the eight equations reduces naturally 
to a system of six first order equations as the displacement 
component U and the corresponding force T. vanish due to the 


symmetry of the applied load. 
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CHAPTER 3 


STIFFNESS ANALYSIS 


Sole Introduction 


This chapter describes how the standard methods of 
structural analysis can be employed to obtain a complete solution 
in terms of the stresses and the displacements at any point within 
a complex structure composed of many elements of any type of 
shell of revolution, using the eight fundamental equations derived 


in Chapter 2. 


3.2 Influence Coefficients 


If a perfectly elastic structural component, supported 
against rigid body motion, is acted upon by a set of forces F), 
Bote Ooo aie es FO See POINeS fe Le Wont soe end 1b Cine ade On CO 
this set of forces, intermediate loads are applied simultaneously 
to the structural member, the induced deformation d. at the point 
i due to the forces ES = 1, 2, 3, «.«5 nm) can be expressed in 


the form 

25 O 
d = i a Pad Seok 
where 7 bees is the additional deformation at point i due to the 

i 


intermediate applied loads and ae represents the deflection d. 


due to a unit. value of es The elastic constants ate are 
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independent of the magnitude of the applied forces ei Eq.: sels 
which represent the linear relation of deformations and forces, 


can be written in matrix form as 


i eee (Ath) +. {4} ne 


where the square matrix [A] is known as the flexibility matrix. 


For the same problem the relation given in Eq. 3.2 can 
also be written in a form which represents the equilibrium 


conditions for the structure as 


{F} [K]{a - d°} 


Il 


[K]{d} + {F°} 323 


where the square matrix [K] is defined the element stiffness 
matrix. The element ee is defined as the force ae due to unit 

oe fe) fe) fe) 
displacement at d.. {F } is the column vector ia 5 oes, ee ee ot 
which are the equivalent forces that replace the intermediate 
applied loads (equivalent in the sense that the work done during 
any incremental deformation approximates the work done by the 
actual applied load). These forces are numerically equal to 


the so-called fixed end forces. 


As stated above, while establishing the flexibility 
matrix the structure is assumed to be supported against rigid 
body motion, a condition not necessary for the stiffness matrix. 


For the latter case the structure can be free to move as a rigid 
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body when a set of nodal displacements is applied. The stiffness 
matrix, thus obtained, is called "the direct stiffness matrix" 
which is singular. The elements in any column represent a force 


system in equilibrium. 


The well known advantage of the direct stiffness method 
is that for an assembly of structural elements, the total 
structural stiffness matrix can be easily formed by superposition 
of the individual stiffness matrices irrespective of the boundary 
conditions. The boundary conditions are considered only after 
assembly in the actual solution of the system of equations. 

This permits the consideration of different boundary conditions 


while the total structure stiffness matrix remains unaltered. 


3.3 Shell Element Influence Coefficients 


By a procedure similar to that used for a beam element, 
the stiffness matrix for an element of shell of revolution can be 
obtained. Since the shell element is a rotationally symmetric 
element, the nodal points are replaced with nodal circles (Fig. 
3.1). The forces and displacements can be expressed as the 


amplitude of a harmonic number along the nodal circles. 


Consider the shell element shown in Fig. 3.1. The 
eight fundamental stress resultants per unit length at the two 
edges of the element are shown acting in the positive sense 
according to the shell theory sign convention. Should the edge 


s = a of the element be subjected to a unit value of displacement 


42. 


sn oar ‘ 
Oy > 


saeeeatae oft a heriqre/etess: 


‘ie y 
iS¥em peana thas: P9Stch cles 


astot & JROSSIVEST Hao 


it 


S9dtS fits SG22 silanehe revi? ita 
taton ‘oils aan Binehts nanan th 1 si a8 ay 

a 

ott tacriaane f Heirs G2 Yiress Ad ids aie orsitanpar de 
nsctt 4dt De Set oaGe ote stil tue < paiaiteae | 
%5 1 00 ore 4h) sit p- Sah eiraras ie cya wt Re 


_ 


yolk 
Om 1) pe dae eas v3 abl aveie 
aS 
1 Ho! VEE ba eed 913 2 tere nition heii, ants 
era uw eerthmer xiv jen Beenie sos Rabo 


oe 


Lay 
Pra 


danidiok 12BO9) £ sti il 


: 26) } 7 


» etuena ls meen. s sod (beeu\ IER mei ileal i: ia 


r 
by ey 


70 Heo WoL ditove®. Io Liane Re naa ul ai tage 8 

Shas sunmvs, WbLon igeapy oe Wi welo Pee end son 

a | arty 2esiasts febao shite ee Alga. Si Ala! satan 0 
ki grit ef, pueoaitgxe 58 gayi 4 sriadaetga’ & ‘te, 870% 1g ig a 
seSloat4 Uabiog 3" Hi ‘Od Fs seq ae ie ns | 


a ae | 
rane ad] , 1... etd, ei 
L i] 


2 | “ed sis 3 aang Jha x49 sanding aesate, ot ra 
5 i aS 


x ser 
if | Ben ¥ Y Peer oe shept ais md be 7h, oct on tei’ me id : 7 7 
| " 5 r , lll epi eesvas2 gic ern, vie 
anon nid 5 Re bnpsediie “28 
wed mf 
a ‘ : ad a : * ft i an . a y 
: 7 ha, Lae wa vA 9 ; ee | y ior 
i i ¢ : y - 7 


acy it ee ‘ ve 


in the direction of W while preventing any other displacement 
of the two ends of the element, the forces per unit length that 
would be required to maintain equilibrium of the element are 
the influence coefficients for a unit displacement Wo: If one 
writes the degrees of freedom of the element in the order Bo Le 
ia oe B, W 


b V, UL» the element stiffness equation (Eq. 3.3) 


can be written as 


x 
M B M° 
sa a sa 
S W oe 
sa a sa 
N V n° 
sa a sa 
T U ‘a 
sa a sa 
= [K] + O 3.4 
M 
Na BD sb 
O 
ay ahs : sb 
O 
N 
eo i sb 
fe) 
Te 
hep nn sb 
fe) fe) fe) fe) : 
Whee Me sacy is Ne SL 2 are the fixed end stress resultants, 
Ss Ss Ss 


at the ends a or b as subscripted, due to the applied loads or 


thermal gradients. 


In order to establish the stiffness matrix and the 
fixed end stresses, the predetermination of the relationship 
between the edge displacements and the edge forces is necessary. 
The stiffness coefficients and fixed end forces for a beam element 
can be evaluated by the well known method, such as consistent 


deformation, or solving the pertinent differential equations. 
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If an analytical solution of the relevant differential equations 
for a specific geometry of shells is known, both influence 
coefficients and fixed end forces due to the applied loads can 
be found. Such a solution is not available for an arbitrary 
element of shells of revolution with variable thickness and 
subjected to arbitrary loadings. However, the forces can be 


evaluated by solving the basic differential equations numerically. 


Sat, -SOLUtiOn, of the Governing System of Equations 


For any given shell element with given geometry and 
applied load, the governing system of equations (Eq. 2.24) is 


in the form 


ere): = [Ay }ty(,y} + {Bey} B15 


where eye is a vector of the eight dependent variables, 
four displacements and four corresponding forces, at a 
particular location s = constant; aye is a vector of the 
derivatives of the eight variables with respect to the 
eoordinate s.i.e., d/dsty.y}5 [A(.)] is the coefficient 
matrix relating the variables and their derivatives and 
consists only of functions of the shell thickness and the 
geometrical parameters of the location s; and, {By ey} is 

a vector of the inhomogeneous terms in the equations and is a 


function of the applied loads and the location s. 
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The géneral solution..of Eq. 3.5. consists of two 


parts: 


He) The solution of the homogeneous part of the equations, 
i.e., the differential equations when all loading 
terms, {Bey}, are set equal to zero. This homo- 
geneous solution involves the evaluation of eight 
constants of integration as the result of eight 
boundary conditions at the two discontinuous edges 
of the shell. 

2) The particular solution of the equations in which 
all loading terms are considered. This solution 


does not depend upon the boundary conditions. 


Therefore, the general solution that satisfies the 
governing system of equations together with the appropriate 
boundary conditions at the two edges of the shell element can be 


written as 


ee = {fh 


where {P y} represents the particular solution that satisfies 
Ss 


( 


the equation. 


} + {B( } car 
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and, the.y} is the homogeneous solution that satisfies the 


equation. 


a wh ae [A(,y] they} 3.8 
Consider now the solution of Eq. 3.8 for an element 
which spans the region b > s > a. Since Eq. 3.8 represents eight 
first order linear differential equations, the solution should 
contain eight arbitrary constants of integration. Let these 
arbitrary constants of integration be the eight (arbitrary) 
boundary values which can be imposed at edge "a'', and denote 


these values by {c}. Then at edge "a" 
ieee we et 3.9 
bubStiatuting into Eq. 3.3 yields, for s =a 


{ } e [Ay] {c} cake 


2 (s)'s = a s=a 


Integrating this numerically, as an initial value problen, 


allows the value of he at any point s > a to be determined as 


= [H Lic, cea) 
the.) [ (s)! 
where (HY | represents the matrix arising from the integration of 
Ss 
the [A ] matrix along the length of the element. 


(s) 
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Since the matrix [Hy] is independent of the initial 
aonditions ({c}), it is a property of the element and may be 
interpreted as follows. As {c} is arbitrary, assign a unit 
value to one component, say component oe and set the others 
equal to zero. Then the jth column of [Hy] represents the 
values of the.y} for a unit value of the jth boundary condition 


at Ue 


It is apparent that when s = a, Eq. 3.11 must reduce to 
Eq. 3.9 in order to match the arbitrary boundary conditions, and 


enc Gol -cclumist tbe thesidentdtyematrix,, ie. 


(a) 


[HY .y! = [I] a1 


Eq. 3.12 may be considered to be a "boundary condition" on the 


numerical integration of the matrix [H]. 


Turming mow to the solution of Eq. 3./, the equation at 


S = a may be written as 


{ } a ie teks st {Brey} Seats. 


es s=a 


where {c*} represents an arbitrary set of initial values of 


{P }. Numerical integration of this equation yields 


(a) 


{P(.y} = [H(, ]te*t + {Q(.y} coed 


where [ ] is the matrix that arises in Eq. 3.11 and {Q.y} 


ve 


is a vector arising from the integration of the inhomogeneous 


terms. Since the particular solution is any solution which 
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satisfies the inhomogeneous equations it is adequate to 


select 
{c*} = 0 
in which case Eq. 3.14 reduces to 


{Py} = {Q,y3 


and consequently the general solution, Eq. 3.6, becomes 


ty (gy} = [Hy ]icd + {Qgy/ 


For s = b, Eq. 3.17 becomes 


1b)! = [Hp Ity cay} + {Qe} 


in which each column vector of [H] represents the variables at 


"b" corresponding to each unit variable applied at "a" 


absence of any applied loads. The vector {Q} represents the 


She a 


Dako 


variables at '"b'' corresponding to zero displacements and stresses 


resultants at "a'’ in the presence of the applied loads. 


3.5 Shell Element Stiffness Matrix 


The column vector ty (gy represents the vector of 
Eq. 2.24 consisting of the four displacements as defined by 


Eq. 2.15.1 and the four stress resultants as defined by Eq. 
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Let 


a 
ty (ay) = : B19 
2) 
and 
(D 
Ria.) = . 3.20 
(b) i, 


Eq. 3.18 can be written in partitioned form as 


where 
D represents the four displacement variables, 
F represents the four stress resultant variables, and 
Qa> Qe are the displacement and the stress resultant 


parts of the particular solution respectively. 


The total matrix appearing in Eq. 3.21 is usually referred to 
aona.. Cranerer matrix’) 123). 78a. 3.2) can be: expanded to-form 


two equations as follow 
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solving Eq. 3.22¥for the.vector < dD, F and substituting 


PLO Bde. 3.23, yields 


F D 0 
a a 
=e ek a eae ta + Saas 
"bl Pipa ted lei Nee 
or 
F D ae 
a a a 
= [K] ses 3225 
ES Ps oes 


It can be seen that each column of [K] represents the stress 
resultants at each end for a unit displacement applied at one 
end while the other displacements are restrained and {F°} 
represents the stress resultants corresponding to the totally 


restrained boundaries. 
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3.6 Stiffness Matrix Sign Convention 


In the derivation of the element stiffness matrix and 
the fixed end stresses, the sign convention used corresponds 
to that generally used in shell theory as given in Fig. 3.1. 
As a result, the stiffness matrix will have some negative 
elements on the main diagonal. This can be corrected by adapting 
the so called "stiffness matrix sign convention”. This sign 
convention is shown in Fig. 3.2. It can be seen that the positive 
direction of the top normal in plane force No» the top tangential 
shearing force T o> the bottom moment M. and the bottom transverse 
shear S have been changed to the opposite direction. Therefore, 
the stiffness matrix and fixed end stresses in Eq. 3.25 are 


to be premultiplied by the diagonal matrix 
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3.7 Stress Resultants and Displacements at Intermediate Points 


For a shell structure composed of a number of elements, 
the element stiffness matrices and fixed end forces are evaluated 
and assembled in the master stiffness equations of the structure. 
Boundary conditions are imposed, and the displacements at the 
boundaries of each element are obtained from the solution of 
the master equations. By substituting the final known boundary 
displacements of each element into the corresponding element 
stiffness equation (Eq. 3.25), the primary stress resultants at 
the element boundaries can be obtained. Thus the correct boundary 
conditions, displacements and primary stress resultants, for 
each element are known. In order to evaluate the displacements 
and the stress resultants at any desired number of intermediate 
points within the element, the correct boundary conditions at 
one end of the element, say the end at s = a, are used as initial 
conditions in integrating the governing set of equations (Eq. 3.5). 
At each intermediate point, the secondary stress resultants 
(which were eliminated from the governing equations) can be 


evaluated, first by evaluating the derivatives of the displace- 


ments using Eq. 2.18. 
iD} = PALM DE [A2]{F} + {A3} S226 
and then substituting into Eq. 2.20 


{Fa} = [cl]{p°} + [c2]{pD} + {c3} cae 
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A simple check on the results of the integration is 


Se, 


that the displacements and the primary stress resultants at the 


termination end of the element should agree exactly with the 


known boundary conditions at this end. 


3.8 Transformation from Local to Global Coordinates 


along the generator, are presented in the direction tangent to 


Displacements and stress resultants, at any point 


the meridian at this point and the direction perpendicular to it. 


Due to possible discontinuity of the meridian curve at a junction 


between two elements, it is necessary to transform the influence 


coefficients and the fixed end forces at this junction, to a new 


coordinates system. 
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It is simplest to adopt the direction of 


According to the stiffness matrix sign 


convention, the transformation equations may be written as follow 
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{D} and {F} represent the displacement and forces respectively 
and the subscripts L, G represent local and global coordinates. 
o. is the angle measured from the global axis of rotation to the 


meridian at the point s = i (Fig. 3.3). 


3.9 Modification of Stiffness Coefficients to Account for 
Eccentricity 

In many cases the middle surfaces of two elements 
which meet at a node do not coincide at the same point (Fig. 3.4). 
A transformation of the stiffness coefficients and fixed end 
stresses to a common reference point is then necessary before 
assembling the element stiffness matrices into the master 
stiffness matrix. Eqs. A.10 relate the displacement components 
of a point at a distance 42 from the middle surface to the 
displacement components of a point on the middle surface lying 
in the same plane. By expanding Eqs. A.10, by means of Fourier 


series, one can write the following equation 
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where Eq. A.7.1 has been used to eliminate the derivative of W 


with respect to s and the subscript z refers to the point at a 
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distance z from the middle surface. Inverting Eq. 3.30 results 


in 
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From the work equivalence requirements, the relation between the 


stresses at the two points can be written as 


oe = [EC]*{F} 3°33 


Eq. 3.32 and 3.33 are used to transform the displacements and 
stresses of a point on the middle surface to a point at a distance 


Zz from the middle surface. 
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CHAPTER 4 


EXAMPLE APPLICATIONS 


fo) watroduct.ion 


A computer program, named SASHELL, has been developed 
to perform the stiffness analysis of segmented shell structures 
based on the theory presented in the preceding chapters. The 
logic flow of SASHELL is outlined in Sect. 4.2 and listing of the 
program is included in Appendix C. The results of the analysis 
of two example applications, using SASHELL, are presented in 


this chapter. 


The first example is the pinched cylinder. The exact 
analytical solution for a long cylinder pinched by a symmetrical 
Piveumrercicial line Lond (ric. 452) Ls known [273° pp. 4725 lo. 
pp. 280]. Finite element solutions, using (48 x 48) element 
stiffness matrices [11] and (24 x 24) element stiffness 
matrices [5,7], were obtained for the case of a cylindrical shell 
loaded by diametrically opposed concentrated loads (Fig. 4.2). 
This solution was compared [2] against an analytical solution 
based on the inextensional deformation theory (i.e. neglecting 


entirely the strain in the middle surface of the shell) [27, pp. 
501-506] . 


The second example is the analysis of hyperboloid 


natural draft cooling tower under the action of wind load and 
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its own weight. Finite element analysis, using the computer 
program SORIII [14] and using conical shell elements to approxi- 


mate the geometry [24] are known. 


ie Logic Flow of SASHELL 


In this section, the organization of the computer program 
SASHELL, which can serve as a summary for the solution technique 
of segmented shell structure as described in the preceding 
chapters, is outlined. Details of the required input are given 
in Appendix B. 
1) The structure is divided into elements, each of which 
is a simple type of shell of revolution, and which 
are connected ‘along nodal circles or "nodes". A 
concentrated load applied at a point along the meridional 
coordinate must be treated as a load acting on a node 
connecting two elements. 
2) Nodal coordinates, system connectivity information 
and element types, properties and loading conditions 
are determined. The problem control parameters are 


established and input (Subroutine READIN and LOADIN). 


The program now performs the following operations: 

1) Each element is examined. If a coefficient which 
depends on the element geometric parameters (see the 
limitations in Sect. 5.3) exceeds a certain limit, 
the element is divided into subelements (segments) each 


of which satisfies this limit (subroutine SEGMNT). 
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The connectivity of the system is altered to include 

the new intermediate nodal points. Geometric parameters 
and loadings values at each segment boundary are 
calculated (subroutine SEGEOM). 

The number of the structure's degrees of freedom are 
established. 

If the external applied load on the structure is 
symmetric, the program discardsstep 5 and goes to 

step 6. 

If the load is non-axisymmetric, the total number of 
points in the meridional and circumferential direction 
in the structure is determined and the "results" array 
is initialized. 

The harmonic number, n, is set equal to zero. 

The structure stiffness matrix and load vector are 
initialized. The nodal load coefficients, if any, 

which correspond to the harmonic n are added to the 

load vector. 

The stiffness analysis starts by calculating the stiffness 
matrix and fixed end forces for each segment in the 
structure (subroutine STIFAN). 

The geometric parameters and external applied loadings 
at the desired number of integration points in the 
segment under consideration are calculated. Dead weight 
of the segment, if required and if n = 0, is superimposed 


(subroutine PLSEG and DLSEG). 
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UE) 


2) 


13) 


14) 


15) 


16) 


The initial conditions at the starting edge of the 
Seement, as'stated in Sect. 3.4. (Eqs. 3.12 and/3.15), 
are set. The governing equations (subroutine FLUGGE) 
are integrated, using a fourth order Runge-Kutta 
method, over the desired number of points to obtain 
the transfer matrix of the segment (subroutine 
RNGKT). 

The segment stiffness matrix and fixed end forces are 
evaluated from the transfer matrix obtained in step 10, 
as deseribed. in. Sect., 325, Gsubroutine.STIFIX). 

The results of step 11 are saved (subroutine STORE1). 
The segment stiffness matrix and fixed end forces are 
modified to correspond to the stiffness matrix sign 
Convention, as, statedi in Sect. 3.6. 

The stiffness influence coefficients and fixed end 
forces are modified to account for nodal eccentricity, 


if any, as mentioned in Sect. 3.9. (subroutine ECCNTR). 


Gri, 


If required, the segment stiffness influence coefficients 


and fixed end forces are transformed, due to discontinuity 


of the meridian at the node, to the structural global 
coordinates, as stated in Sect. 3.8 (subroutine 
GLTRAN) . 

The segment stiffness influence coefficients are 
assembled, with respect to the structure's degrees 

of freedom, into the master stiffness matrix. The 


segment fixed end forces are subtracted from the 
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corresponding values in the load vector (subroutine 
SLORE Y's 

Steps 9 to 16 are repeated for each segment in the 
Structure. 

The boundary conditions are imposed on the master 
stiffness equation (subroutine BOUNDC). 

Segments edge displacements is obtained, using 
Gaussian elimination algorithm to solve the master 
stiffness equation (subroutine SOLVER). 

The stiffness analysis, for harmonic number n, is 
completed. The displacements and stress resultants 

at the desired number of intermediate points in each 
segment are to be evaluated (subroutine SRADSP). 

The segment edge displacements, obtained in step 19, 
are transformed, if necessary, to the segment local 
coordinates and to account for nodal eccentricity 
(subroutine GLTRAN and ECCNTR). 

The known segment edge displacements are substituted 
in the corresponding segment stiffness equation, saved 
in step 12. The primary stress resultants at the 
segment boundaries are evaluated, as described in Sect. 
3.7 (subroutine STORE1). 

The shell equations are integrated, using the initial 
conditions calculated in step 22, in order to determine 
the displacements and primary and secondary stress 
resultants at the intermediate points within the segment, 


as described in Sect. 3.7 (subroutine RESULT). 
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24) If the load is symmetric (i.e., the required number of 
harmonics is zero), the results of step 23 are printed 
out and the program goes to step 26. 

25) If the load is non-axisymmetric (i.e., the required 
number of harmonic is greater than zero), the displace- 
ments and stress resultants are calculated at the desired 
number of points in the circumferential direction. The 
results are superimposed in the "results" array 
(subroutine STORE2). 

26) Steps 21 to 25 are repeated for each segment in the 
SELUCTUTe. 

27) The harmonic number, n, is increased by one. Steps 
7 to 26 are repeated until n is equal to the required 
number of harmonics. 

28) The results saved in step 25 are printed out and the 


program stops. 


4.3 Pinched Cylinder 


The governing differential equation for a circular 


cylindrical element is 


d2 d“wW ae of 7 
FED (KS = ts iW P. 


where x measures the distance along the axis of the cylinder 
with the origin (x = 0) at midlength of the element. When the 


thickness of the shell is constant, this equation reduces to 


64. 


aA 
sv 


) 


vedabn ba stwpar wll cate io 


y) 7h Co Gewese eae na 9 a 


Lupas 


> 
< 
3°30 
uf | 

of 
Le ® 


jabantie Pr eee aren gantt 


ne alee 


fel ae: ey Sas Bis) ae 


ebuparr Gilt ay) yia Peery ere ae bak se 
add iors carte epee Oh Cimersad, So 7 


nie - Sot idiappe ites lanl 


ak quire gt Oye ‘aexgone ca pra 


) las 


det4siene oh acon tah? g8Q132 sing a 


4 ~ ta, 


inthe dee Meh bahay ‘te : a ; 


oe Miah anette ne istenalueun oye me tuenty iT a. 

eaasnnte, jonhaukrisiuay 9 ai 

+ yt sabiasy Tot borkedghaetm Cfi02 09 ane " 
; seen nora. 

rv eh Bie bornes eae acti anes wit, aa 

Sx) i i) pay is Tor dag bern areal ae is : 

Le 


Hy Ris TG, ied st Sie Co ie 
sey 


an . M ri a nt 


) 


a) fi 


ae 


‘ i : . Mises aus 
} - g / : 
a! >; ; ; 7 
uy \ a \p T : ; on + sigs “ir 
a 4 
= ih i) , 
| : i wi 2 


65. 


d*w 4 Py 
sTeop cae W = a Lo? 


Poevnach ene flexural risidity, K, is 


hi Be" 
Sit CL aes) 4.3 


ena the coefficient A is defined as 


ba ieodel ar Ve} 
For the special case of a long cylindrical shell pinched by a 
live load P uniformly distributed along a circular section 
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The stress resultants may be expressed in terms of derivatives 


of the displacement W in the form 
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Numerical values of these expressions, at several points along a 
cylindrical element 20 feet long, 8 feet in diameter and with 
constant wall thickness equals to 1.24 inches, are shown in 
Peolesg- | torea live load P = 1 kip/fr: (Outputvof the analysis 
of the same element using SASHELL is included in Appendix D. 

The results are presented graphically in Fig. 4.3. The results 


of the two solutions are identical up to the number of significant 


figures contained in the output. 


For the case of a cylindrical shell pinched by two 
concentrated loads as shown in Fig. 4.2 the concentrated load is 
approximated as a live load (see Sect. 5.4) that has a value 
of 1.0 kip/ft at the loaded points and zero at a short distance 
from the load points (Fig. 4.4.1). The loading function is 
expanded in Fourier coefficients and the results of the analysis 


of each harmonic are superimposed. 


The results of this loading case are affected by 
the number of harmonics considered in the analysis, which will be 


diccuseed insect. 5.4. 


If the circumference of the cylinder shown in Fig. 4.4 
is divided to 36 intervals and the load is approximated as shown, 


the equivalent concentrated load is 


0.698 kips 
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TABLE 4.1 Stresses and Radial Displacement in a Pinched 
Cylinder as per Eqs. 4.6 
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The analytical solution for a concentrated load, based on 


inextensional deformation theory [27, pg. 506], is of the form 
z eee 4.8 
(ns) , 


where L is one-half the length of the cylinder. Eq. 4.8 yields 


(W) Ss POsos2 x LON ULL 


a -3 
(Woionsg = 70-351 x 10° ft 


A comparison of these values with SASHELL values shown in Fig. 
4.4.1 indicates good agreement for 8 = 1/2 but considerable 
discrepancy for 9 = 0. The SASHELL solution predicts a displace- 


ment 17% larger than the inextensional solution at 6 = 0. 


Now referring to Fig. 4.4.2, the cylindrical shell 
example as choosen by Cantin and Clough [7], is shown. For the 


same number of intervals (36), the equivalent concentrated load is 


P = 0.0721 kips 


The deflections for this case, from Eq. 4.8, are 


= O)xola se 10> fe 
(W) 5 0.614 x 


Sik itp oa 
(W) gon 9 0.565 x 
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The ratios of the SASHELL displacements are 1.10 and 1.14 to 
those of the inextensional theory for 0 = 1/2 and 9 = 0, 


respectively. 


4.4 Hyperboloid Natural Draft Cooling Tower 


Large capacity power plants generate a substantial 
amount of operational heat that requires dissipation. One of the 
major structures in these power plants is the natural draft 
cooling tower in the form of the shell of revolution. The 
tower utilizes its height to create the necessary air flow in 
order to cool a large volume of water in a minimum land area. 
Hyperboloid cooling towers are the most preferable shape [14] 
when compared to conical or cylindrical shapes from the aero- 


dynamic point of view. 


ary ee Geometry of the Tower 


The middle surface of a hyperboloid shell is shown in 
Fig. 1.7. The surface of a hyperboloid of revolution may be 
classified as a non-developable surface, which means that the 
surface will not tend to flatten out under load. This surface 
is generated by the rotation of a hyperbola about a vertical 


axis. The geometrical equation can be written as 


: sinh ie =} 4.9 


ce 
ine) 
o 


in which r is the horizontal radius, x is the vertical coordinate 
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measured from the origin at the throat of the shell, a is the 
throat radius at x = 0, and b is a constant in which the ratio 


b/a equals the slope of the asymptotes to the hyperbola. 


The principal radii of curvature are given by 
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The expression for the derivative of r, with respect to the 

coordinate s (see Eq. 2.5.1) can be obtained by differentiating 

Eq. 4.10.1 "to yield 
-3 
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The middle surface of the hyperboloid tower considered 
in the following analysis is shown in Fig. 4.5. The shell is 
355 ft high and is supported by columns evenly spaced on a 
circular base of 290 ft diameter. The throat of the tower is 
165 ft in diameter and is located 60 ft below the top of the 


shell. The thickness varies from 30 inches at the bottom level 
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of the shell to 6 inches at 25 feet elevation from the bottom. 

In the top 10 feet of the shell, the thickness also varies from 

6 inches to 24 inches. Other than the top and the bottom regions, 
the shell thickness remains constant at 6 inches. The increased 
thickness at the top provides a stiffening effect that reduces 
radial deformation under wind load [24]. The bottom ring at 


the base acts as an equivalent deep beam bridging between columns. 


The constant b of Eqs. 4.10 can be calculated by 
substituting the values for r = 145 at x = 295 and a = 82.5 
mieound. 4.9, which yields b = 204.1 ft. For the purpose’ of 
comparing the results with References 14 and 24, the following 


concrete properties are used 


Young's Modulus Beta 4 x 10° pet 
Poisson's Ratio yeaa 0. 5 
Specific Weight Agen USC Lb fhe 


4.4.2 Load Description 


The tower is analyzed for gravity load and wind 


pressure load. 


The gravity load, which is symmetric with respect to 


the coordinate 9, is determined from the following equations 
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Pe = q sin 6 Gk e 3 


where q is the intensity of weight per square area of the surface. 


The wind load, which is the governing factor in the 
design of the cooling towers [14], is based on the ACI-ASCE 
Committee 334 recommendations [1]. The following equivalent 


static pressure distribution is used 


q(H,9) 7 G Co K 930 eB 


where the factors in this equation are described as follows: 


a) q(H,8) is the equivalent static normal pressure on the 
surface of the tower at a location defined by coordinates 
H,9. H is the vertical distance measured from the ground 
level, 8 is the circumferential angle measured from the 
windward meridian. 

b) G is the dynamic gust factor which accounts for the 
overstress due to the tower response to various time 
variations in the wind pressure. 

e) Cy is the coefficient of wind pressure distribution 
in the circumferential direction. 

d) Ky is the exposure factor which establishes the 
vertical profile of wind pressure, which in turn, 


depends on wind speed and roughness of terrain. 


i.e., K. can be evaluated from the equation 
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where H is the height measured from the ground level; 
Hg is the gradient height above which the wind velocity 
is assumed constant and ranges from 900 ft for open 
country to @500.¢t for center of Jarge cities: ca issa 
constant, depending on the terrain roughness, and 
ranges from 1/7 for open country to 1/3 for center of 
harce cities. 

e) q30 is the basic wind pressure (psf) and is equal to 
the dynamic pressure of the free stream of wind at 
30 ft above the ground level at a given site. It may 


be computed from 


qao = 0.00256 V3, 4.14 


where V39 = wind velocity (m.p.h.) at 30 ft above the 


ground level. 


The ACI-ASCE Committee 334 has suggested the normalized wind 
pressure distribution in the circumferential direction as shown 
in Fig. 4.6.1... .As in Reference 24, assuming Hg = 900 ft, q = 1/7 
and V39 = 100 m.p.h., the variation of the wind pressure with the 


height can be taken as shown in Fig. 4.6.2. 


4.4.3 Analysis Conclusions 


The results of the analysis of the hyperboloid tower, 


loaded as described above, are represented graphically for 
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dead load in Figs. 4.7.1 to 4.7.3, and for wind load in Figs. 
4.8.1 to 4.8.10. Since the actual condition at the bottom of 
the shell is only partially fixed [24], the results of the 
analysis for two boundary conditions, fixed and hinged at the 
bottom, are shown in these figures. Excellent agreement is 
observed between the displacements and stress resultants shown 


and the results of References 14 and 24. 


The approximation of the geometry with a series of 
cones [24] does not affect the results for this particular 
geometry of the tower. However, the meridian curvature R, varies 
only from -0.002 f£t7! at the throat to -0.0003 ft~! at the bottom 
of the hyperboloid shell, and is equal to zero for a conical 
element. For shells with larger R, curvature the conical 


segment approximation may not be as accurate. 


4.5 Comments on Results 


The examples of this chapter have been selected to 
test the ability of SASHELL to analyze different types of shells. 
The pinched cylinder solutions are common test problems because 
of the difficulty of achieving solutions for concentrated loads. 
The hyperboloid cooling tower is an illustration of a shell with 
negative Gaussian curvature under complex loading conditions. 

It may be concluded that SASHELL is capable of yielding good 


results on a wide variety of shell problems. 
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Fig. 4.1 


Fig. 4.2 


x in eq. 4.3 


“Pinched Cylinder’’ Circumferential Line Load 


“Pinched Cylinder’” Two Concentrated Loads 
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——— Sashell Solution 
- Analytical Solution [27] 


. 
a Kel Ol ratte 


0.1249 k ft/ft 


0.5 k/ft 


far 7 Mg 
—_—_— 0.0366 k ft/ft 


Ng 
Ses 3.995 k/ft 


Fig. 4.3 
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CHAPTER 5 


LIMITATIONS 


boi Wntroduction 


Three factors affect the solution technique presented 


in this thesis. These are summarized as follows. 


i) Singularity of the governing equations at the apex. 
2) Stability of the numerical integration process. 
5) Convergence of Fourier expansions. 


Each factor is discussed separately in the following sections. 


5.2 Singularity of the Governing Equations at the Apex 


itethesshell+has a pole (i.e., r = 0), coefficients 
in the governing sets of equations (Eqs. 2.14, 2.18 and 2.20) 
become singular. This is consistent with classical shell theory. 
A simple way to handle this situation is to choose the boundary, 
s = 0, not at the pole, but a very short distance away and 


then impose the boundary conditions at s = 0 as follows [6]. 


1) For harmonic number n = O 


Ds) For harmonic number n = 1 
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3) For harmonic number n > lL 
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DPS Stability of the Numerical Integration Process 


Runge-Kutta fourth order integration method is used, 
in SASHELL, to integrate the basic set of equations (Eq. 2.24). 
This method is very well known and has the advantage of being 
self starting i.e., it needs only the information available 


at the preceding point. The method can be interpreted as follows: 


1) The derivative is evaluated at the starting point of 
the interval. 

2) The above derivative is used to obtain an approximate 
ordinate to determine an approximate derivative for 
the midpoint of the interval. 

s) The above derivative is used to obtain a second 
approximation of the derivative at the midpoint of 
the interval. 

4) The above derivative is used to obtain an approximate 
ordinate to determine an approximate derivative for 
the end point of the interval. 

5) A weighted average of the above four derivatives is 


taken to determine a total increment in the function 


for the whole interval. 


Analytically, this can be defined as follows: 
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in which £(Y,, h,) is the value of the function at point i, and 


Ah is the step size of the interval. 


It has been found [25] that the direct integration 
methods, when applied to a shell problem, suffer a complete 
loss of accuracy when the generator of the shell exceeds a critical 
length. The reason for this phenomenon is explained clearly 
in Reference 18. The general solution of the governing equations, 
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Because of the exponentially decaying behavior of the stresses 
and displacements, it is observed that the coefficients of [H] 
increase in magnitude in such a way that if the length of the 


shell element is increased by any factor n, then these coefficients 
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increase in magnitude, approximately, exponentially with n. 

For example, consider an element spans the region a SoS et De 
For some prescribed edge conditions at "a", we expect the 
corresponding solution at '"b'’ to become smaller and smaller 
when the element ab is increased in length. (i.e., ie. small, 
Ob) large and eS has a prescribed value). The longer the 
element, the larger [Hoy] and the smaller ty pyt. The only 


] 


way to get a small value for ty py} is for the elements of [Ee 
to subtract out and that is the reason’ that at some critical 


length of the element all significant digits of [Hopy ls in 


Eq. 5.3, are lost and so is the accuracy. 


The loss of accuracy cannot be avcided by choosing a 
fine mesh for the integration. By taking more steps "partial 
instability" arises [8]. This means that the numerical solution 


deviates from the actual solution as we take more steps. 


In Reference 18, the critical meridian length is 
limited with a length factor AL < 3 - 5, where 
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in which 
L is the length of the meridian of the shell, R is the 


minimum radius of curvature, and h is the thickness of the shell. 


In the author's opinion, this is very conservative 
limitation. The loss of accuracy, using the computer program 


SASHELL, does not arise until the length factor AL exceeds 25. 
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To demonstrate this, Table 5.1 shows some results of the analysis 
of clamped pressurized cylindrical shells, as obtained by the 
computer program SASHELL, with various lengths and different 
number of steps of integration. From the symmetry of the problem, 
the end moments and the absolute values of the end shears should 
Peweduat. lt is obvious, from Table 5.1; that® the solution 1s 
not affected much by the number of steps adopted for the inte- 


gration and it begins to break when 


Ie) AL exceeds 20 and the number of steps is less than 20. 


2) AL approaches 29 for any number of steps. 


As a conservative limitation, the upper bound of AL shall 
not exceed 25 and the number of points of integration (NP), needed 


for convergence, is limited as follows: 
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This limitation does not affect the efficiency of 
this technique. It can be observed from Table 5.1, that the 
computation time (CPU) required for the analysis is proportional 


to the number of steps adopted. 


5.4 Convergence of Fourier Expansions 


The load, when varying with the circumferential 


coordinates is expanded in a Fourier series. The analysis is 
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carried out to determine the stresses and displacements ever ywhere 
within the structure. When the load is represented "exactly", 


the solution converges. 


Theoretically, the number of harmonics required for 
an arbitrary periodic function to be represented exactly, by 
Means (of..Fourier series. is infinity, i.e., 


ee) oe) 
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n=0 n=l 

However, the load can always be described at a sufficient number 

of points to satisfy our engineering judgement of representing 

the actual loading conditions. Therefore, our concern is to 

examine the convergence of a function known only at a set of 

discrete points. If these points are equally spaced, say 2N 


points, taken over the interval 0 < Oo < 27,. the spacing ts, 
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where the coefficients are to be determined in such a way that 


the integrated squared error over the interval of length 27 is 


least, then the requirement is 


96; 


i ~ ‘y pe 
0 bt i 
» oe " : mares ; ao 
; n ee min) ney 


Hwy-07 |S Sap. BTS shal a 


, yey) Relea ae ae eiuk eae ie 


= ; ii’ 3 hi iL eso. ag anes nas boat 
7 ; : nay 

(Sa toy Fateh peep mte ipod hana be 
hi Yer 200. are. eee yee) anatdtboos wokberd) fy) 


Vio wéioas smal 3o0o7 @, ty sbiio ase 


a 


i 
is 


Vii Pee ihe s 3 oe Panton seats, —— wut 


nat sage oll ae h- > y Lavmasas. etd 300 et GD e 
: ' 0 ; : A a 7 3 


Te 


"e 14 ae 


sn at myo ati 
4 | X Ngee 
ii it Le 9 alee ha 


As. . Nowe, nh 


24 Che gigas nee a 


oF". 


277 


M M 
| Fg) mp: AL cos! nl; Ss Rarsiin eno = min 
0 n=0 n=L iS 


Since the function F Og) is described only at the points defined 
by o. (Eq. 5.6), we have 2N independent data points which are 
sufficient to determine the coefficients of 2N terms of an 
approximation in the form of Eq. 5./. When M < N, Eq. 5.8 can 


be rewritten as follows 


2N-1 M M 2 
. fig on LAY Cos? noe 2), San no. = min 
n=0 ae n=0 2 , n=l 


where f isthe value of the function £ a pcheypodn bas 


(84) (03) 
The solution is obtained when the partial derivatives of the 
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tovzero (15, pe. 446-457]... The coefficients are im the torm 
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The calculation of each coefficient is independent of the 
calculation of the others and is independent of M as long as 
Ms N. When M=N, the least-squares criterion (Eq. 5.9) 
becomes equivalent to the requirement that the "best" 


approximation of the function (Eq. 5.7) is obtained [15]. 


The loading cases of the examples discussed in Chapter 4 
ere shown in Fig. 5.1 and.5.2. The effect of the. number of 
harmonics on the loading associated with these problems is 
considered in the following. Fourier coefficients of the wind 
pressure load on the hyperboloid tower are included in Table 5.2.1. 
Since the load is symmetric with respect to 8 = 0, the sine 
coefficients vanish. The results of superimposing the 12 
Poerticients, using. Eq... 5./,\.are: shown in Table 5.2.2. Lt 
can be seen that these values approximate the load function to 
a very close agreement. Now examining the coefficients shown 
in Table 5.2.1. The 9 and the subsequent coefficients are 
small in comparison with the other coefficients. The magnitude 
of the gth coefficient is in the same order of the 10th with 
opposite sign and so for the 11th and 12th, Therefore we expect 
that the approximation will be reasonably accurate if the series 


is terminated at the gth harmonic. The results of an 8 harmonic 


approximation are shown in Table Sie ees 


Fourier coefficients of the two concentrated loads of 
the pinched cylinder problem discussed in Sect. 4.3 are shown 
in Table 5.3.1. It can be seen that the odd cosine coefficients 


as well as the sine coefficients vanish for this case of symmetry 
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with respect to 9 = 0 and 9 = 7/2. The approximation with 18 
harmonics, for the loading function described at 36 points, 
represent the load more accurately when compared with the 12 
harmonic approximation (Fig. 5.2 and Table 5.3.2). Therefore 
one can approximate this load by superimposing the results of 


TE, 1s iso ets: LO. 


As a general conclusion, one can consider N harmonics 
in the expansion of a load described at 2N points. A termination 
of the higher harmonics or exemption of a harmonic number in 
the series can be decided upon by examining the coefficient of 


each loading function as a special case. 
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8 Load Value for Load Value tor 
12 Harmonics 8 Harmonics 
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TABLE 5.2.2 Effect of Number of 
Harmonics on Representing 
the Wind Pressure Load 
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Fourier Coefficients 
of Two Diametrically 
Opposed Concentrated 
Loads 
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TABLE 5.3.2 Effect of Number of 
Harmonics on Representing 
the Two Concentrated Loads 
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Fig. 9.2 Fourier Approximation for Two Diametrically Opposed Concentrated Loads 
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CHAPTER 6 


SUMMARY AND CONCLUSIONS 


In this study, a theory has been generalized. A 
computer program has been developed for the elastic analysis of 
axisymmetric segmented shell structure of general geometric 
configuration. General arbitrary loadings have been considered. 
Applications for a number of loading cases and elements have 


been presented. 


The reduction of the governing partial differential 
equations of the classical shell theory to a set of eight first 
order ordinary differential equations involves only straight forward 
algebraic manipulations. The classical theory chosen as a 
foundation for the theory presented herein, Flugge's theory, is 
considered one of the most accurate theories available in the 
literature. The approximation in formulating the basic governing 
equations of this theory are such that the theory may be considered 
to be exact. The adaptation of the general form of these basic 
equations allows the geometry of any type of shell element to 
be considered, without any theoretical approximations which may 
be suitable for specific dimensions and not for others. 

Variation of the shell thickness along the meridian can be 


accounted for with accuracy comparable to that for constant 


thickness. 
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The loads, when approximated by means of Fourier series, 


may be represented with sufficient accuracy by including a number 


of harmonics, within practical limits. This approximation 
allows the consideration of general arbitrary types of loadings 


in a simple manner. 


The excellent agreement between the results of the 
problems presented in this thesis and the results of other known 
solutions, analytical or numerical, demonstrate the accuracy of 
the solution technique and the reliability of the numerical 


integration process used. 


It may be concluded that the solution technique 


presented in this study is simple and general for the elastic 


analysis of shells of revolution. The accuracy is consistent with 


a good shell theory. The geometric limitations are of no 


importance from the practical point of view. 


Further development, by applying this method to the 


study of free vibration and elastic buckling criteria of shells 


of revolution is possible. 
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APPENDIX A 


SHELL THEORY 


This appendix presents the analytical development of 


the basic equations, used in this study, as given by Flugge [10]. 


A.1 Geometry of Shells 


The geometry of a shell is defined by specifying the 
form of the middle surface and the thickness of the shell at 
each point. The surface of a shell of revolution is generated by 
the rotation of a plane curve about an axis in its plane. This 
generating curve is called the 'meridian". The intersection of 
the surface with planes perpendicular to the axis of revolution 
are "parallel circles". Two coordinates s, 8 are required to 


describe any point onthe middle surface of the shell: 


a) S measures the distance to the point along the meridian 
from the intersection of the middle surface with the 
axis of) rotation, or from a datum parallel vcircte: 

b) 6 is the angular "distance" of the point from a 


datum generator. 
A third coordinate z is required to measure the distance along 
a normal to the middle surface. 


The radii of curvature of a shell of revolution are: 
a) r is the radius of curvature of parallel circles 


b) r,; is the radius of curvature of meridian 
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r2 is the length of the normal between any point on the 


middle surface and the axis of revolution. 


The following fundamental geometrical relations can clearly be 


seen in Fig. A.l 
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epresents the angular distance of the point under 


tiontf£romithesaxis of ‘rotation. 


From Eqs. A.l one can write 


“x = cosd tee ara h 
CL te eh A as iad A oe A.2.2 
ds Y. 


For simplicity the derivatives with respect to s and 
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ee Equations of Equilibrium 


If a shell element is cut out by two meridians and two 
parallel circles, each pair infinitely close as seen in Fig. A.2, 
the element is stressed by ten stress resultant components which 
must be in equilibrium with the external applied load. These 


stress resultants are: 


No» Ny = normal in-plane forces per meridional and circum- 
ferential unit length, respectively. 

Nig? No, = in-plane shear forces per meridional and circum- 
ferential unit length, respectively. 

Qo» Qy = transverse shear forces per meridional and 

circumferential unit lengths, respectively. 

Mo> My = meridional and circumferential moments per unit 
length, respectively. 

Mog? Mp = circumferential and meridional twisting moments per 


unit length, respectively. 


The sign convention for the forces is shown in Fig. A.2 


Referring to the three orthogonal axis S, 0, z, one can 


obtain six equations of equilibrium 
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» M, is the sum of the moments about the i axis 


(i = S5 he, z) 


The six equilibrium equations in terms of the ten stress resultants 
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equations of equilibrium are not enough to determine the ten 
stress resultants. The shell element is four times internally 
statically indeterminate and therefore, an analysis of the 
deformation of the shell is required in order to obtain the 


solution. 


Ao Strain-Displacement Relations 


The displacement vector of a point lying on the middle 
surface of an element may be described by its three orthogonal 


components, defined as: 


W = the displacement component in the radial direction, positive 
when it points away from the centre of curvature. 

V = the displacement component in the meridional direction, 
positive in the direction of increasing the coordinate s. 

U = the displacement component in the direction of the tangent 
to the parallel circle, positive in the direction of increasing 


the coordinate 0. 


In addition, the auxiliary variable 8, which represents the angle 
by which an element of the meridian rotates during deformation 


may be expressed in terms of the displacement components, Fig. A.3, 
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Referring to Fig. A.4, one can obtain the relationship between 
the strains and the displacement components of a point on the 
middle surface. 

The meridiofial strain is 

2 _ Elongation of the line element ds 

Ss ds 
ery. A.8.1 
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The hoop strain is 


_ Elongation of the line element rd0 
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The shear strain, which is the change of the right angle between 


the two-line elements: ds and rd0 (Fig. A.4.3), is-equivalent to 
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Since the displacement is assumed to be very small in 
comparison with the principal radii of curvature, all products of 


two displacement components have been dropped in deriving the 


foregoing equations. 
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One may use the preceding equations to express the 
relationship between the strains and the displacement components 
of an arbitrary point at a distance z from the middle surface 
by simply replacing W, V, U with the displacement components 
Of this point Wo vo ee and replacing the radii ri, reo with 
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Therefore, the strains at a distance z from the middle surface 
in terms of the displacements W, V, U at the middle surface can 


be obtained by substituting Eqs. A.10 into Eqs. A.9 to obtain 
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Hooke's law relates the strains to the corresponding 


within the elastic limit. 


measured from arbitrary level. 


index notation, as 
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stresses in linearized form, as long as the stresses remain 
If T is the change in temperature 


Hooke's law may be written, in 
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Oo, as the normal stress in the i diréction 


i 
Ti is the shearing stress in the i plane and 
j direction 
E 
G = OL +0) A.13 


The modulus of elasticity, E, Poisson's ratio, v, the coefficient 
of thermal expansion, a, and thus the shear modulus, G, are the 


material constants. 


As in the theory of plates, except in the immediate 
Wieinity Of concentrated forces, the stresses in z direction are 
small in comparison with the stresses in s, 9 directions and 
their influence in Hooke's law may be neglected, i.e., it is 


assumed that 
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Therefore, the stress-strain relations can be written as 
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wo Elastic Law 
The internal stress resultants can be determined by 
integrating the stresses through the shell thickness (Fig. A.6). 
They are defined as 
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It should be noted that the expressions for Cer Qe have been 
omitted as they are equal to the integral of the shearing 

stresses in z direction. Also the minus sign in these expressions 
correspond to the positive directions assumed for the stress 


resultants as shown in Fig. A.6. 

The expressions for 05> Oy» T X26 (Eqse- #215) Gan be 
entered into Eqs. A.16. Then Eqs. A.11 are substituted for the 
displacements and the integration with respect to z is performed. 
The results in the stress-resultant-displacement relationships. 


ye wy . U' + V.cos@ + W sind 
1 


ZA 
I 

oO 
=| 


POSES vat sag AA ee a bp +2 
ry 1 ies 


ae Wi cue eee Bs S) re Aelia t 


[ u'_ + V cosh + W sind 4 Vy. +H 
1 


K 5 ey Eo al Sat a2 (as) cosd =e W sind 
ar Bo igy, 


Vecibetes W° cose | - (ior) oe D ae Arad 
1G 


{ - ae Ge es : 
bt ae 
; ' ‘v5 (oe 
i. ” 
4 % ri ‘2 ” ‘ 
4 mt Ve Sy = 
* ; 
ir or A j 
cz 
rh 8 € aU 
u >» @ 
ay 
: amt rie % ar AGE at 
7 


i} 4 le | butoge - x 


a 


im. nae oe 
mr 


.) 


sO 


Os 


D ea : Vy [e ue v' - : =| 


K buat fs = ry | 6 pbk = 165) 
+ — oS ES, U —S—_____*. 


Vicosh 2%, - T1 4 \yyee pent en ee 


aiteny) 1a) Tl 
Vv ae oe 

v $y — 5 | 1 +) oO KK Ti, 
igi oa 


ae. 


Are fens 


Aiea 


A.1/.6 


—— - 


Rhos ae, ti Ae id oe a 


a | f . 9 


4 7 at na 


2 


Lie ae 2w' U° 


i 2 


and the flexural rigidity, K, is defined as 
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If a linear variation of the temperature T through 
the thickness is assumed, Eqs. A.19 can be integrated by parts to 


yield 
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in which T° and qt are the temperature measured at the outer 

and inner face of the shell respectively. If the term t/ry, is 
neglected when compared with unity, the subscript k disappears 

from Eqs. A.20 and as T1 can be defined as the average temperature 
measured on the middle surface and the temperature gradient 


respectively. 


By substituting Eq. A.6 into Eqs. A.17 and changing 


the subscript s to $, Eqs. A.1/ reduce to the elastic law in 


Reference 10, pg. 322. 


Eqs. A.5 and A.l/ are the governing equations for a 
shell of revolution. The sixth equation of equilibrium (Eq. A.5.6) 
i= identically satisfied 1f rm = re. dfheretore; the five 
remaining equations of equilibrium (Eqs. A.5) and the eight 
equations of the elastic law (Eqs. A.17) are 13 equations in 13 
unknowns (the three displacement components and ten stress 
resultants). Theoretically, one can solve for the stresses 


and displacements at any point in the shell using these equations. 
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Fig. A.1 Meridian and Parallel Circle of a Shell of Revolution 
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Fig. A:3.1 Meridional Rotation Due to Displacement V 


Fig. 4.3.2 Meridional Rotation Due to Displacement W 
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Fig. A.4.1 Meridian of a Shell Before and After Deformation 
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Fig. A.4.2 Parallel Circle Before and After Deformation 
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Fig. A.4.3 Change of the Right Angle Between Line Elements After Deformation 
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Fig. A.6 Stresses Acting on a Shell Element 


APPENDIX B 


USER'S MANUAL FOR PROGRAM SASHELL 


Program SASHELL computes stress resultants and displace- 
ments for axisymmetric branched segmented shell structures due to 
their own weight, external applied loads and differential 


temperature variation (along the meridian or circumference). 


In the present stage of development, the program is 
capable of analysing five types of shells of revolution of variable 
thickness. These are cylinders, circular plates, spheres, cones 
and hyperboloids of revolution. Loadings may be symmetric or 
non-axisymmetric with respect to the collatitude coordinate and 


may vary along the meridian. 


The analysis procedure is based on the theory presented 
in this thesis and the program logic flow outlined in Sect’ 4.2. 
A complete listing of the program is given in Appendix C. Input 
of the problems discussed in Chapter 4 and output of the pinched 


cylinder for a concentrated line load are given in Appendix D. 


The input to SASHELL consists of several types of 


input cards. Certain card types may be repeated as required. 


A typical explanation of a card type consists of the 
card type, a descriptive name indicating the nature of the data 
being entered and the format for the data on that card. This is 
followed by a symbolic line of input which, in turn, is followed 


by definitions of the input variables. 
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Limitations on SASHELL, due to the dimensions of the 


arrays in the program, are outlined following the explanation 


of the input cards. 


TYPE Ls 


EYRE 23 


IPRINT 


TITLE CARD (Format 10A8) 


80 


AN IDENTIFIER STRING 


One card which contains any title for the problem 


ANALYSIS CONTROL CARD (Format 414,F7.0) 


4 8 12 16 Z3 


IPRINT NPCR LDC BETA 


Print control parameter 

If IPRINT = 0, the output will contain an echo check 
of the completed data. The loadings will be expanded 
in Fourier series, if required for the analysis, and 
the coefficients will be printed out. 

If IPRINT = 1, the output will contain the echo check 
of the input data and the final results. 

If IPRINT = 2, the output will contain the echo check 
of the input data, the results of the analysis for 
each harmonic and the superimposed final results. 

If IPRINT = 3, full output including intermediate 


values will be printed out. (Used for checking 


purposes only.) 
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NPCR 
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BETA 
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Number of points along the element meridian for 
which the Runge-Kutta integration process is used. 
NP should not be specified less than 21 (see Limitations). 
Number of the circumferential points at which the 
final results are required. 

If the loads are symmetric or antisymmetric with 
respect to the meridian passing through 6 = 0, NPCR 
is the number of points along half the circumference 
(0, 7). If the loads vary randomly in the circum- 
ferential direction, NPCR is the number of points 
along the full circumference of the element (0, 27). 
If the loads are constants in the circumferential 
direction (symmetric), NPCR is equal to one. 

Dead load control parameter. 

If LDC = 0, dead load is excluded from the analysis. 


1, dead load is evaluated and superimposed 


free 
on the external applied loadings. 
Maximum element length coefficient. BETA should 


not exceed 25 (see Limitations). 


STRUCTURE DATA CARD (Format 213,4F12.0) 


3 6 18 30 42 54 


Number of elements 


Number of junctions between elements (nodes). 
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EG : Global modulus of elasticity. 

PUG : Global Poisson's ratio. 

GAMG : Global specific weight. 

TKG : Global coefficient of thermal expansion 


NOTE: If the structure consists of elements of different 
materials, the global properties are to be omitted and the 


structural data card specifies the number of elements and nodes 


only. 
TYPE 4: NODAL DATA CARDS (format 14,2F10.0,414) 


One card is required for each node. 


4 14 24 28 a2 36 40 


L : Node number 

XCOOR(I) : Giobal X coordinate of node I along the axis of 
revolution directed downward from the top of the 
Sisbucmthice:. 

RCOOR(I) : Radius of the parallel circle passing through node I. 
RCOOR(I) should not be specified as zero (see 
Limitations). 

IDF(I,J) : Identification of the jth degree of freedom at node I. 
J = 1,4 for the rotation of the meridian (8), the 
radial displacement component (W), the meridional 
displacement component (V) and the circumferential 


displacement component (U), respectively. When 
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IDF(I,J) = 0, the corresponding degree of freedom 
is not restrained. 
When IDF(I,J) = 1, the corresponding degree of 


freedom is restrained. 


TYPE 5: ELEMENT DATA CARDS 


Two cards are required for each element. 


First card Format (516,4F10.0) 

6 12 18 24 30 40 50 60 70 
12 : Element number 
rt 1) : Element type 


If IT(I) = 1, element I is a cylinder 


2, welementel, LS 2 COne OF. a 


Lee LL) 


circular plate. 


If IT(I1) = 3, element I is a sphere of which rj = ro. 
Tf IT(1) = 4, element 1 1s a sphete of which ci to. 
If IT(1) = 5, element I is a hyperboloid of revolution. 

RGtd 1) : Node number at the top of element I. 

BGCh Se) : Node number at the bottom of element I. 

NIP (1) : Integer to indicate the number of intermediate points 


in the element at which the final results are not 
required. The number of equally spaced points at 


which the final results will be printed out are 
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NP-1 
NIP (T) + 1. A value of NP-1 will be assigned to 


NIP(I) when input is zero. 

PCL.) : Element thickness at the top. 

TCE. ) : Element thickness at the bottom. 

ECCE,<L) : Eccentricity of the top node from the middle 
surface of the element at the top. 

ECC 2) > Eccentricity of the bottom node from the middle 


surface of the element at the bottom. 


NOTE: The eccentricity is defined such that the radius of the 


parallel circle passes through the middle surface of the element is 


eee one pe eee 


Thus, EC(I,J) is positive when directed inward from the node to 


the middle surface of the element. 


Second card (Format /F10 20) 


10 20 30 40 50 60 70 


HPCN(1,1) : Radius of curvature of the meridian for spherical 
element of type 4 (i.e., ri # r2), or throat radius 
of a hyperboloid element (type 5). For elements 
other than of type 4 and type 5, HPCN(1,I) is 0.0. 

HPCN(2,I) : Angle in degrees measured from the axis of revolution 
to the top edge of the spherical element of type 4, 


3. WEPONC IT) 
or hyperboloid constant in which the ratio HPCN(1, 1) 
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equals to the slope of the asymptotes of the 
hyperbola. For elements other than of type 4 and 5, 
HPONG225) as Od. 

HPCH(3,I1) : Angle in degrees measured from the axis of revolution 
to the top edge of the spherical element of type 4, 
or global X coordinate of the throat of the hyper- 
boloid element. HPCH(3,1) is equal to 0.0 for 
elements other than of type 4 and 5. 

GAMA (1) : Element specific weight. Specified if different 
from the global specific weight, otherwise GAMA(I) = 0.0 
or blank. 

ECT) : Element modulus of elasticity. Specified if 
different from the global modulus of elasticity, 


otherwise E(1) = 0.0 or blank. 


PU (IL) : Element Poisson's ratio. Specified if different 
from the global Poisson's ratio, otherwise PU(I) = 0.0 
or blank. 

TCOEF(I) : Element coefficient of thermal expansion. Specified 


if different from the global coefficient otherwise 
TCOEF(I) = 0.0 or blank. 
NOTE: For a segmented structure which does not include spherical 
or hyperboloid elements and for which the other elements of which 


the structure consists have the same material properties, this 


ecard is a blank card. 
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TYPE 6: LOADING SPECIFICATION CARD (Format (615)) 
5 10 D5 20 ZD 30 
NEL NJL NHL NTL NHPL NHIN 
NEL : Number of externally loaded elements in the structure. 
NJL : Number of loaded nodes in the structure. 
NHL : Maximum number of harmonics required for the analysis 


including the zero harmonic (i.e., NHL = maximum 
harmonic number +1). 

NTL ;- Loading type character. 
If NTL = 0, loading is symmetric and the analysis 
is required for the zero harmonic only. 
If NTL = 1, loading is non-axisymmetric in the 
circumferential direction, input is provided at a 
number of discrete points along the circumference of 
the shell, and the analysis is required for the cosine 
coefficients of Fourier series only (i.e., loading 
is symmetric with respect to a meridian passes 
through 9 = constant). 
If NTL = 2, loading is non-axisymmetric in the 
circumferential direction, input is provided at a 
number of discrete points along the circumference of 
the shell, and the analysis is required for both are 
cosine and sine coefficients of Fourier series. 


If NTL = 3, loading is nonaxisymmetric along the 
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along the circumferential direction and input 
directly as cosine coefficients of Fourier expansion 
only. 

NHPL : Number of points along the circumference of the 
shell at which the loads values are described. 


These points are defined such that the circumferential 


: at ee : 
coordinate is NHPL wherey ints. Ohi clase Wie. eeeNEeL. — 1. 
NHIN : Integer which defines the increment in the harmonics, 


starting from the zero harmonic, to be specified 
when the analysis is required for a harmonic number 
COSWNHIN, 2 NHINoweis ese NHL ~ oL)\ » The: programscets 
NHIN = 1 when it is specified as zero or blank. 


NOTE: If the load is described at NHPL points, NHL can be 


NHPL 
Zz 


IPRINT = 0, the Fourier coefficients of the input load are 


specified as and NHIN can be set equal to one. When 


obtained with the echo check of the input data. Then, the user 
may decide, upon examining these coefficients, on the final values 


of NHL and NHIN (see Chapter 5, Sect. 5.2.3). 


TYPE 7: ELEMENT LOADING CONTROL INPUT CARDS (Format 615) 


One card for each loaded element 


5 10 is 20 25 30 
£1, Ct) TROAD. (ELC) 78) bh eee 


i : Integer takes the value of 1 to NEL 
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Number of the loaded element. 


ILOAD (LL (I) ,K) 


NOTE: 


Identifier for loading type on the element LL(T) 

in the order: 

ILOAD(*,1) for loading in the direction tangent 

to the meridian (s). 

ILOAD(*,2) for loading in the direction tangent 
tothe paralbel Circle (6) . 

ILOAD(*,3) for loading in the direction perpendicular 


to the tangent to the meridian (z). 
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ILOAD(*,4) for temperature at the shell exterior face. 


ILOAD(*,5) for temperature at the shell interior face. 


0, no load of type K is applied. 


If ILOAD(*,K) 
If ILOAD(*,K) = 1, the applied load of type K is 
constant along the meridian and to be specified at 
one end of the element only. 

If ILOAD(*,K) = 2, the applied load of type K varies 
linearly along the, meridian.and,is to, be specified 
at the two ends of the element. 

If ILOAD(*,K) = 3, the applied load of type K 


varies as a second degree function along the meridian 


and is to be specified at the two ends of the element. 


Element loads are input in the element local coordinates. 
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TYPE 8: ELEMENT LOADING CARDS 


This type of card depends upon NTL and is classified 
as TYPE 8-A, TYPE 8-B and TYPE 8-C for NTL equal “to0.et pore. 


and 3 respectively. 


TYPE 8-A: ELEMENT SYMMETRICAL LOADING CARDS (Format 2F10.0) 


This type is required if NTL is equal to zero. Number 
of cards required, for each loaded element, is equal to, and 
input in the order consistent with, each non-zero term in the 


corresponding row in ILOAD array. 


10 20 
ACEL(K1) ACEL (K2) 
ACEL(K1) : Magnitude of the load at the element top. 
ACEL(K2) : Magnitude of the load at the element bottom. 
TYPE 8-B ELEMENT ASYMMETRIC TABULATED LOADING CARDS (Format 8F10.0) 


This type is required if NIL is equal to 1 or 2. 
Number of cards required, for each loaded element, is equal to 


NHPL for, and input in the order consistent with, each non-zero 


8 


term in the corresponding row in ILOAD array. 
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Paik 
NHPL 


(i = 0, 1, 2, ..., NHPL-1) along the circumference 


Magnitude of the load, at the points defined by 


of the element, and for which Fourier expansion is 


required. 


TYPE 8-C ELEMENT LOADING FOURIER COEFFICIENTS CARDS 


(Format 8F10.0) 


This type is required if NTL is equal to 4. Number of 


NHL 


cards required, for each loaded element, is equal to 8 


input in the order consistent with, each non-zero term in the 


corresponding row in ILOAD array. 


AL (K) : Cosine coefficients of Fourier expansion for 
Harmonics OF b,  2ilyss.ce NHR 
TYPE 9: NODAL LOADING CONTROL INPUT CARDS (Format I5) 


One card for each loaded node (NJL cards) 


5 


ss 4 
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ak - Number of loaded node. 
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TYPE 10: NODAL LOADING CARDS 


This type of card depends upon NTL and is classified as TYPE 
10-A, TYPE 10-B and TYPE 10-C for NTL edualstou0, pecore-aand ios 


respectively. 


TYPE 10-A: NODAL SYMMETRICAL LOADING CARDS (Format (4F10.0)) 


This type is required if NIL is equal to zero. One card is 
required for each loaded node. 


10 20 30 40 


! if 


! 
ANGLCIE, 1). MET) =) 1,4 


ANA Giy 1): Magnitude of the nodal load, in the global structure 
coordinates, in the following order. 
II = 1, meridional couple, Mo: positive when rotating in 
the clockwise direction. 
II = 2, force normal to the tangent to meridian at the 
corresponding node, Ss positive in the outward direction 
from the axis of revolution. 
II = 3, force in the direction tangent to the meridian, 
No» positive when directed downward parallel to the axis 
of revolution. 
II = 4, force in the direction tangent to the parallel 
CLTEC Le, Lae positive when directed in the anticlockwise 
rotating direction around the structure. 

NOTE: The sign convention of the nodal loads as defined above is equiv- 


alent to the stiffness matrix sign convention as described in Sect. 3.7. 
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TYPE 10-B: 


NODAL ASYMMETRICAL TABULATED LOADING CARDS 


This type is required if NTL is equal to 1 or 2 and 


consists of the following cards, for each loaded node. 


(1) 


(2) 


Asymmetric nodal load control input card. Required 
to identify which to the nodal loads, classified as in 
Type 10-A cards (M_, So No> T.)> is applied an to be 


input. One card (Format 414) is required 


4 8 16, 420 


1 t 
KEDCK)) , Ko = 1,4 | 


K : Integer which takes a value of 1 to 4 and 
represents the nodal forces, as defined in 
TYPE 10-A cards, in the order Mo> S > N. and 
T o> respectively. 

KEDCK) va! ih equal tto.zero, nodabrfencetoh the type K 


is not applied < 


Ef sequal ctoels-nodaliforcevanh thestype K is 


applied. 
Nodal loading magnitude card (Format 8F10.0) 
NHPL 
Number of cards required is equal to 3 for, and 


input in the order consistent with, each non-zero term 


in, ILD 


10. 20 .s30). 405550 60270080 
' 
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1 t ! t ' 


WCLY 4, P= Ll, NEEL 
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TYPE 10-C: 


Magnitude of the nodal load, at the points defined 


2771 


by NEPL (i= 05/4,02, 23.5 )0HPIA =) alone the 


circumference of the nodal circle, and for which 


Fourier expansion is required. 


NODAL LOADING FOURIER COEFFICIENTS CARDS 


This type is required if NTL is equal to 3 and NJL is 


greater than 0. It consists of the following cards, for each 


loaded node: 


(1) Asymmetric nodal load control input card of the same 


type described in TYPE 10-B. 


(2) Nodal loading magnitude card (Format 8F10.0) 


NHL 


Number of cards required is equal to ae LOL, and 


input in the order consistent with, each non-zero 


Lerm in KED; 


AL(N) 


Cosine coefficients of FOURIER expansion for 


harmonics 0, 1, 2, .-- NHL-1l. 
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LIMITATIONS 


Number of elements NE 

Number of nodes NJ 

Number of integration points NP 

Number of circumferential points NPCR 

Number of harmonics NHL 

Number of circumferential points at which loadings 
is specified NHPL 

Full band width NHB 

Total number of Fourier coefficients in the problem 
for elements loading, associated with either cosine 
or sine factor and defined by NHL x NEL 

Total number of Fourier coefficients in the problem for 
nodal loadings, associated with both cosine and sine 
factor and defined by 8 x NHL x NJL 

Number of points in a non-axisymmetrically loaded 
structure at which the final results are required, 


and which can be calculated as 


(NP-1) 
——+~ + ] 
) NS; x (yp. j 
i=l 8 
where NS. = number of segments into which element i is 
a 


ip — 
subdivided and can be obtained by BETA 4/ Wd = fee 


20 


aah 


51 


11 


CN ENO NG Ne NG 


20 


# 40 


# 200 


? 200 


7 200 


in which, L, is the element length, r, is the element radius, t, 


is the element thickness and v is Poisson's ratio. 
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The above limitations are due to the dimension 
statements in the program SASHELL as listed in Appendix C. The 
dimension of the corresponding arrays may be modified to change 


these limitations. 


Limitations on the theory employed in SASHELL are as 


follows (see Chapter 5). 


1) RCOOR should not be specified as zero. 

es) The length factor coefficient BETA should not be 
specified greater than 25. 

S)) The number of points of integration NP should not be 
less than 21 when BETA is specified less than 20, and 
not less than 31 when BETA is greater than 20 and less 
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